The Global Positioning System (GPS) has become very valuable for communication, navigation, and for geolocation. The GPS system includes about twenty-four spacecraft in 12-hour orbits at about 11000 nautical miles above the Earth. Each spacecraft includes a precision clock and a signal transmitter. A GPS receiver receives signals from several GPS spacecraft and measures the time of emission of the signal from each spacecraft. It then uses this information, together with knowledge of the location of the GPS satellites derived from the GPS signal, to calculate its own geolocation and time. Navigation systems exploit this knowledge to control vehicles and to target weapons. The GPS system also supports the synchronization of network systems such as communication systems.
One important aspect of GPS performance is the ability to predict the GPS SV (space vehicle) position and clock state, as the broadcast message is based upon prediction. Currently the prediction interval is about 24 hours. The dominant error in the range prediction error is the instability of the SV's clock. By enhancing the stability of the SV clock, range prediction errors are reduced and the accuracy of the GPS system is improved.
A clock includes a periodic process and a counter. A simple clock produces a number in a register, which number corresponds to how many counts or how many periodic events have taken place since the clock was initialized. Most physical processes are not perfectly regular in time due to small physical noise processes that vary over time. Hence, clocks will deviate from “perfect” time, and the difference between absolute time and time measured by a clock will behave as a random process and the variation of the error will grow with time.
Clock stabilities have been characterized classically by Allan variances. To measure the Allan variances at specific time samples (Δt), a sequence of phase errors (xk) are measured (the difference between phase readings between two clocks, one being a more stable clock standard). This is used to compute a sequence of double differences of phase errors (ak=xk+2−2xk+1+xk). Then the Allan variance, avar, is computed using the following formula:
                              a          ⁢                                          ⁢                      var            ⁡                          (                              Δ                ⁢                                                                  ⁢                t                            )                                      =                              1                          2              ⁢              Δ              ⁢                                                          ⁢                              t                2                            ⁢              n                                ⁢                                    ∑                              k                =                1                            n                        ⁢                          a              k              2                                                          (        1        )            
To compute the Allan variance curve, the value of the Allan variance is computed for a whole range of values of (Δt). Typical Allan variance curves for three clocks are given in FIG. 2, plotted in log-log space. One should notice there are three regions of the curve with different characteristics. For shorter sample times the curve has a negative 1 slope; for long sample times the curve has a +1 slope. For intermediate sample times, the curve is flat. Each of these regions is due to a different noise source: phase random walk, flicker, and frequency random walk, respectively. The stability characteristics are related to the ability to predict clock states. The following relationship holds between prediction errors due to these error sources and the Allan variance curve:P(Δt)=a var(Δt)Δt2  (2)
Thus it can be seen the ability to have low clock prediction errors is related to low stability characteristics.
To exploit clocks in an application, such as GPS, it is necessary to generate state space representation of clocks. The minimal state space includes two states: phase and frequency. In certain clocks the center frequency changes linearly with time. For those clocks, a third state, drift (the derivative of frequency) may be modeled. As part of the modeling process, statistics of process noise are modeled to match the stability characteristics. For example, the model:
                              [                                                                      x                  .                                                                                                      y                  .                                                              ]                =                                            [                                                                    0                                                        1                                                                                        0                                                        0                                                              ]                        ⁡                          [                                                                    x                                                                                        y                                                              ]                                +                      [                                                            u                                                                              v                                                      ]                                              (        3        )            Represents a clock with phase random walk u and frequency random walk v (x is the phase error; y is the frequency error). By choosing the covariances of u and v properly, the model will match the Allan deviation curves for those noise sources. There is no simple model that will generate flicker noise; such models have developed and usually involve more states (e.g. Van Dierendonck developed a model with a third state). The most general two state model is given by process noise with the covariance:
                    Q        =                  [                                                                      q                  11                                                                              q                  12                                                                                                      q                  21                                                                              q                  22                                                              ]                                    (        4        )            Q is a positive definite 2×2 matrix. Using the empirical definition of “Allan variance,” the Allan deviation is given bya var(Δt)=q11/Δt+q22Δt  (5)
Notice there is no term corresponding to flicker noise. If the relationship between prediction errors and Allan variances is used, thena var(Δt)=q11/Δt+q12+q22Δt  (6)This does include a term that looks like flicker noise. This term is not true flicker, as it is not an additional noise source but is due to a correlation between phase and frequency random walks. It can be used to help compensate for flicker noise but is limited by the constraint|q12|≦√{square root over (q11q22)}  (7)due to the positive definiteness of Q.
Prior to 1990, the GPS Master Control Station (GPS MCS) employed the “master clock” for time keeping. In the master clock implementation, a single clock is used to define time, and all other clocks are synchronized to the master clock through tracking these clocks from ground stations. Due to random variations and anomalies in the master clock, the master clock concept had its vulnerabilities. One source of these vulnerabilities of the GPS system is the necessity to keep GPS time close to coordinated universal time (UTC) time, which is Greenwich Mean Time updated with leap seconds. When the difference between the master clock and UTC became too large, then it was necessary to switch the master clock function from one clock to another. Changing the master clock introduced a discontinuity into the GPS service, disrupting the service to those receivers that were in operation at the time of the discontinuity.
In the implementation of the master-clock concept for the GPS, each clock is modeled using parameters in the above models. Since time relative to the master clock is solved for, the relative clock model is solved for. This requires for every clock adding its process noise covariance and the master noise process noise covariance in order to provide the process noise of a single clock. Since all clocks have the same reference, it is necessary to also represent the correlation between the clocks. Thus for example, let Q1 and Q2 represent the process noise for clocks 1 and 2. Let QM represent the process noise for the master clock. The process noise matrix for clocks 1 and 2 relative to the master clock will be given by
                    Q        =                  [                                                                                          Q                    1                                    +                                      Q                    M                                                                                                Q                  M                                                                                                      Q                  M                                                                                                  Q                    2                                    +                                      Q                    M                                                                                ]                                    (        8        )            
In the GPS filter, clocks are only a subset of physical systems that must be estimated. States must be included to model the positions of the GPS orbits, etc. In the standard Kalman filter, two steps are included: a time update and a measurement update. The time update equation is given by: Pk+1=ΦkPkΦkT+Qk  (9)
                                                        x              ->                        _                                k            +            1                          =                              Φ            k                    ⁢                                    x              ->                        k                                              (        10        )            P's represent covariance of that state estimates (x) . Phik (Φk) represents the global transition matrix (Φ for a single clock is given by Equation 37); and Q represents global process noise covariance. These matrices are formed from the individual clocks by forming block diagonal matrices by concatenating the matrices belonging to individual clocks. Equations for these quantities for individual clocks, which constitute a subset of the global quantities, have been presented. Other equations are needed for other states to fill out the global matrices. To complete equations for the GPS Kalman filter specification it is necessary to also specify the measurement processing. Measurements are of time delays between signals being emitted from one source and being received at a receiver. The states of clocks of the source and the emitter appear as states in the observation equation, which have the standard linearized form as:Δzk+1=Hk+1Δxk+1+wk+1  (11)    Δzk+1=Δzk+1−h(xk+1)=predicted residuals    Δxk+1=state space correction (to be computed by the filter)    wk+1=measurement noise.
      H          k      +      1        =                              ∂                      h            ⁡                          (              x              )                                                ∂          x                    ⁢              ❘                  x          =                      x                          k              +              1                                            =          measurement      ⁢                          ⁢      Jacobian      
In the filter, it is assumed that measurement errors are zero-mean white noise with a covariance matrix given by Rk+1. In terms of these quantities, the Kalman filter update is given byΔxk+1=Kk+1Δzk+1  (12)whereKk+1= Pk+1Hk+1T[Hk+1 Pk+1Hk+1T+Rk+1]−1  (13)This correction is added to the predicted state to provide an updated estimate state:
                                          x            ->                                k            +            1                          =                                                            x                ->                            _                                      k              +              1                                +                      Δ            ⁢                                                  ⁢                          x                              k                +                1                                                                        (        14        )            The filter covariance is updated according to the equationsPk+1= Pk+1−Kk+1Hk+1 Pk+1  (15)
The above equations are prior art known as “Kalman filter”. Various versions have been implemented in order to provide improved numerical stability. Bierman has discussed such methods at length (see “Factorization Methods for Discrete Sequential Estimation” (Academic Press, 1977)). It is well known that for the Kalman filter (or factorized methods) to work, it is usually required that all states be observable. In the GPS application (or other clock operations), since all measurements are time differences, there is no way to observe all clocks as independent states. This leads to the introduction of the “master clock” discussed above. In the master clock implementation, corrections to the master clock are not solved for. Clock states for the other clocks are corrections to master clock time.
Since 1990, a composite clock, known as the Implicit Ensemble Mean (IEM), has replaced the master clock arrangement for GPS operations. The IEM composite clock uses a weighted linear combination of clocks, located both in space and on the ground to define “GPS time”. The IEM composite clock decreases the sensitivity of GPS time to the behavior of single clocks. The IEM composite clock applies weights associated with space and frequency estimates provided by a Kalman filter that solves for positions and velocities of the GPS satellites in addition to the clock states of the individual clocks. By changing the weights or contributions of the various clocks to GPS time, the GPS time can be gradually steered to UTC time without significant discontinuities. K. Brown, the developer of the IEM concept of “composite clock”, also points out that composite clocks provide better timing and geolocation accuracy than a master clock. Thus, by using a composite clock, the accuracy with which users determine their locations has improved. The argument for improved timing accuracy is that, if each clock is off by a random error, then averaging a group of clocks produces a more accurate estimation of time. If the covariance matrix of errors in all the clock states (phase and frequency) is available, then the IEM composite clock approach uses maximum likelihood estimation with knowledge of the statistical errors of the estimates of phase and frequency of the individual clocks making up the composite to develop the “best” estimate of time by combining time estimates derived from different sources. This covariance is usually not known, as measurements of clocks with respect to absolute time are not available. Since the covariance is unknown, weights are usually chosen arbitrarily, with larger weights assigned to more stable clocks. In addition to Brown, others have used the notation of averaging a set of clocks to produce a more accurate clock.
In the IEM clock it is assumed that a priori clock estimates {right arrow over (x)}0 and covariance P0 are given. This information is then used to generate maximum likelihood estimates of “composite clock” states as follows. Define the following set of observations{right arrow over (x)}0=Hc{right arrow over (x)}0c  (16)where Hc is a matrix comprised of 1's and 0's, e.g. for a two clock example
                              H          c                =                  [                                                    1                                            0                                            0                                            0                                                                    0                                            0                                            1                                            0                                                                    1                                            0                                            0                                            0                                                                    0                                            0                                            1                                            0                                              ]                                    (        17        )            The maximum likelihood estimate of the composite clock is given by{right arrow over (x)}0c=[HcP0−1HcT]−1P0−1HcT{right arrow over (x)}0  (18)The linear function permits mixing phase states and frequency states, which would occur if there were correlations in the covariance matrix P0. Although these equations, as we indicated above, are applied to the initialization of the filter, it appears new linear combinations could emerge if the composite clock was resolved for at some future time. Brown showed that the linear combination would remain invariant in time under the assumption that the process noise of each clock was proportional to a stability reference standard, a condition that he denoted as “PPN” for proportional process noise. Brown also observed that since all measurements are in time differences, there is now new information that would contribute to an improved estimate of the composite clock. In other words, the composite clock will be driven by a fixed linear combination of outputs of constituent clocks based upon a prior information available when the clocks were initialized
FIG. 1a is a simplified block diagram of a prior art Implicit Composite Mean clock, including individual clocks. In FIG. 1a, composite clock 10 includes first (1), kth (k), and nth (n) clocks illustrated as blocks 12, 14, and 16, respectively, which produce individual clock signals {right arrow over (x)}1, {right arrow over (x)}k, and {right arrow over (x)}n, respectively. Blocks 32, 34, and 36 indicate the weights K1, Kk, and Kn respectively, by which these clock signals are multiplied, and 22, 34, and 26 designate the corresponding multipliers. Multiplier 22 multiplies the output {right arrow over (x)}−1 from first clock 12 by gain K1 from block 32. Multiplier 34 multiplies the output {right arrow over (x)}k from kth clock 14 by gain Kk from block 34. Multiplier 26 multiplies the output {right arrow over (x)}n from the nth clock 16 by gain Kn from block 36. Summing (Σ) Block 38 sums the outputs of these signals to produce the composite clock output {right arrow over (x)}c. In FIG. 1a, the multiplication constants Kl, Kk, and Kn from blocks 32, 34, and 36 are scalar values, selected to provide a minimum variance estimate of the initial clock states, based on extrinsic knowledge. The initial setting of the weights, however, does not explicitly address the issue of the stability of the composite clock.
An alternative embodiment of the composite clock is through the use of a master clock and the weighted sum of phase-lock-loops. The use of phase-lock-loops is a prior art practice in order to provide more accurate phase measurements. A phase-lock loop is a standard circuit that is used in communication and signal tracking to measure the frequency and phase offsets between two periodic signals with nearly the same frequency. FIG. 8 is a high-level diagram of a phase-lock loop 600. Two signals, namely S1 and S2, are multiplied in a multiplier 610, and the product signal is filtered by a low pass filter (LPF) 612. The resulting signal is integrated 614 to provide the phase output of the difference of two input signals.
In the arrangement of FIG. 1b, a master clock 52 controls the electronic sampling or reading of the outputs of the other clocks 52′, 52″. Clock 52′ includes blocks 54 and 58, and clock 52″ includes blocks 56 and 60. In FIG. 1b, the output signal {right arrow over (x)}1 of the master clock 52 is applied directly to the summing (Σ) circuit 38, and is also applied to relative phase counter (RPC) circuits, illustrated as relative phase counter blocks 58 and 60, which control clocks 54 and 56. The output of RPC 58 represents a second clock signal Δ{right arrow over (x)}k, which is applied to a weighting or multiplying circuit 22 for weighting by a factor or coefficient Kk from block 34. The output Δ{right arrow over (x)}n of RPC 60 is applied to a weighting circuit 26 for weighting by a factor Kn from a block 36. The weighted outputs of multipliers 22 and 26 are applied to summing circuit 38 to produce a composite clock signal {right arrow over (x)}c. The master clock 52 measures differences between its time (and/or frequency) and the other clocks 52′, 52″ in the constellation of FIG. 1b through relative phase counters 58 and 60, and produces composite time by adjusting its own time through a correction provided to summing circuit 38 from master clock 52.
In the arrangement of FIG. 1b, the weightings are applied to the relative measurements between the reference clock 52 and the individual clocks 52′, 52″, and the weighted differences are summed with the reference clock signal to produce the composite clock signal. The gains Kk and Kn are chosen generally as described in conjunction with FIG. 1a. If the sum of the weights equals 1, the composite clock time will not be sensitive to the stability of the master clock.
The IEM clock embedded within a filter as described by Brown assumes proportional process noise covariances (PPN). In the PPN case, clocks have similarly shaped stability curves as a function of sample time. FIG. 2 illustrates Allan variances σ2 202, 204, and 206 as a function of sample interval for three individual clocks, such as the clocks of FIG. 1a or 1b. As illustrated in FIG. 2, the values of Allan variance σ2 for plots 202, 204, and 206 include a generally flat portion extending between sample intervals Δt1 and Δt2, together with rising values of σ2 below sample interval Δt1 and above sample interval Δt2. The Allan variance plots 202, 204, and 206 of the three clocks have the same general shape, in that the falling portions, flat portions, and rising portions are generally parallel. Consequently, stability characteristics of every clock can be expressed as a scale factor (λ) times a nominal clock, which serves as a reference stability standard. For example, plot 204 of FIG. 2 can be taken as being the standard (λ=1), in which case the scale factor for plot 202 may be, for example, λ=2 and the scale factor for plot 206 may be, for example, λ=3. In other cases, the PPN assumption does not hold. In particular, PPN does not hold when combining crystal clocks and atomic clocks. Crystal clock and atomic have entirely different stability profiles, e.g. the crystal clock can be more stable for short times, while the atomic clocks are more stable for longer times. Thus, the different portions of the variance curves for non-PPN situations are not necessarily parallel, so proportional scaling cannot be accomplished.
A single clock does not provide an estimate of frequency i.e. time between increments, so that the frequency cannot be known other than through some calibration knowledge of the duration of the clock's “second.” The frequency of a clock is inferred from prior or calibrated knowledge of the clock count. However, if more than one clock is present, the combination of outputs can be combined to update an estimate of relative frequency, based upon relative phase measurements. By using a relative measurement of frequency between two clocks obtained by a phase lock loop process, a better estimate of frequency can be obtained by using the measurements of the relative outputs of the clocks. If the clocks are equally stable, an average of the clock readings will provide a more accurate estimate of frequency than that of a single clock.
When plural clocks are available in a composite clock arrangement, the time determination is achieved by adjusting center frequencies and phases to a common standard through the use of cross polling of the clocks by use of a central clock for reference timing; the reference timing is not necessarily used as part of the composite clock output. With distributed systems, like GPS, a master clock can be used in place of a central polling clock. The difference between the use of a master clock and the cross-polling clock is that the master clock defines “time” whereas the cross-polling clock is used to keep a nominal reference or intermediate time, which is updated by corrections. The reference time is corrected by adding the weighted outputs of the differences between reference time and the other clocks of the composite implicit means. In the GPS system, cross polling is achieved through relative clock measurements between spacecraft clocks and clocks on the ground or by relative clock measurements between space vehicles.
Improved clock weighting for implicit ensemble mean composite clocks is desired that explicitly addresses the need to optimize the stability of the composite clock.